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ABSTRACT 

We study the evolution of supernova (SN) remnants of the first stars, taking proper account of the 
radiative feedback of the progenitor stars on the surroundings. We carry out a series of one-dimensional 
hydrodynamic simulations with radiative cooling, starting from initial configurations that are drawn 
from the results of our earlier radiation hydrodynamic simulations of the first Hn regions. We primarily 
consider explosion energies i?sN = 10 51 — 10 53 erg, appropriate for various types of a single Population 
III supernova, and also explore cases with greater energy to model multiple explosions. In low-mass 
( ^ 10 6 Mq) halos, the stellar radiation significantly reduces the ambient gas density prior to the 
SN explosion. The blastwave quickly propagates over the halo's virial radius, leading to complete 
evacuation of the gas even with the input energy of 10 50 erg. We find that a large fraction of the 
remnant's thermal energy is lost in 10 5 — 10 7 yr by line cooling, whereas, for larger explosion energies, 
the remnant expands even more rapidly with decreasing interior density, and cools predominantly via 
inverse Compton process, adding bulk of the energy to cosmic microwave background photons. In 
higher mass (~ 10 7 Mq) halos, the gas density near the explosion site remains high ( Z, 10 4 cm~ 3 ) 
and the SN shock is heavily confined; the thermal energy of the remnant is quickly radiated away by 
free-free emission, even if the total input energy exceeds the binding energy of halos by two orders of 
magnitude. We show that the efficiency of halo destruction is determined not only by the explosion 
energy but also by the gas density profile, and thus controlled by radiative feedback prior to the 
explosion. We compute the emissivity of the remnants in various energy ranges and examine the 
detectability by future observations. Several implications of our results for the formation of first 
quasars and second-generation stars in the universe are also discussed. 
Subject headings: cosmology: theory — stars:Population III — supernovae 



1. INTRODUCTION 

The cosmic Dark Ages ended when the first stars lit 
up the universe. This 'dawn' of the universe may have 
been accompanied by rather violent events; first super- 
nova (SN) explosions. Various feedback effects on the in- 
tergalactic medium (IGM) are caused by SN explosions, 
and the importance of early SNe can be easily appreci- 
ated by noting that only light elements were produced 
during the nucleosynthesis phase in the early universe. 
Heavier elements such as carbon and oxygen must have 
been expelled by means of SN explosions at some early 
epoch to account for various observations of metal en- 
richment at high redshifts. 

In the standard cold dark matter (CDM) model, the 
first objects are formed in low-mass halos (Couchman & 
Rees 1986; Tegmark et al. 1997; Yoshida et al. 2003). 
Even a single SN explosion can destruct such a 'mini- 
halo' and thus can cause a strong negative feedback ef- 
fect: if the halo gas is completely blown away by a SN 
explosion, the subsequent star-formation is likely to be 
suppressed temporarily in the same region. On the other 
hand, if the shock-heated gas remains trapped within a 
deep gravitational potential, gas cooling and condensa- 
tion may occur even more efficiently thereafter; heavy 
elements processed in the precursor star is dispersed in 
the vicinity of the explosion site, and then the overall gas 



cooling efficiency can be locally enhanced. The strengths 
of these effects and the halo destruction efficiency are 
of considerable cosmological interest, because they may 
largely control the global cosmic star formation rate at 
very high redshifts. 

Recent theoretical studies on the formation of primor- 
dial stars consistently suggest that the first stars are 
rather massive (Abel, Bryan & Norman. 2002; Bromm, 
Coppi & Larson 2002; Omukai & Palla 2003). If the 
first stars are indeed as massive as ~ 200 Mq, they 
end their lives as energetic SNe via the pair-instability 
mechanism (e.g., Barkat, Rakavy, & Sack 1967; Bond, 
Arnett, & Carr 1984; Fryer, Woosley, & Heger 2001; 
Heger & Woosley 2002), releasing a total energy of up to 
~ 10 53 erg. Such energetic explosions in the early uni- 
verse are thought to be violently destructive: they expel 
the ambient gas out of the gravitational potential well of 
small-mass dark matter halos, causing an almost com- 
plete evacuation (Bromm, Yoshida, & Hernquist 2003; 
Wada & Venkatesan 2003). Since the very massive stars 
process a substantial fraction of their mass into heavy 
elements, early SN explosions may provide an efficient 
mechanism to pollute the surrounding IGM (Schneider 
et al. 2002; Yoshida, Bromm & Hernquist 2004; Qian k, 
Wasserburg 2005). 

The physics of astrophysical blastwaves has been ex- 
tensively studied in various contexts (e.g., Ikeuchi 1981; 
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Bertschinger 1985; Vishniac, Ostriker, Bcrtschingcr 
1985; see also a comprehensive review by Ostriker & 
McKee 1988). On a cosmological background, Ikeuchi 
(1981) suggested energetic explosions in the early uni- 
verse as a large-scale star-formation and galaxy for- 
mation mechanism. Carr, Bond & Arnett (1984) and 
Wandel (1985) studied shell fragmentation in a similar 
scenario involving Population III SN explosions in the 
early universe. Dekel & Silk (1986) studied the im- 
pact of star-bursts and associated SN-driven winds on 
the formation of dwarf galaxies, whereas Tegmark, Silk 
& Evrard (1993) suggested a reionization model by SN- 
driven blastwaves. Pre-galactic outflows driven by multi- 
ple SNe are suggested also as an IGM enrichment mech- 
anism at high redshifts (Madau, Ferrara & Rees 2001; 
Mori, Ferrara & Madau 2002). More recently, Bromm 
et al. (2003) carried out three-dimensional (3-D) sim- 
ulations of the first SN explosions in a CDM universe 
and studied the various feedback effects caused by them. 
Wada & Venkatesan (2003) performed grid-based Eule- 
rian hydrodynamic simulations of SN explosions in iso- 
lated small galaxies and showed that the expelled gas 
falls back to the potential well after about the system's 
free-fall time. While these previous works consistently 
showed the destructive aspect of early SN explosions, 
they either simulated a few specific cases or employed 
somewhat simplified initial conditions. The density pro- 
file around the explosion sites can be more complex and 
critically control the efficiency of cooling of SNRs (e.g. 
Terlevich et al. 1992). Therefore, it is important to carry 
out hydrodynamic simulations starting from proper, re- 
alistic initial conditions. 

In our earlier work (Kitayama, Yoshida, Susa & 
Umemura 2004; hereafter Paper I), we studied the 
evolution of cosmological Hn regions around the first 
stars, by solving self-consistently radiative transfer, non- 
equilibrium chemistry, and hydrodynamics. We showed 
that the final gas density profile depends sensitively on 
the halo mass. In 'mini-halos' with mass < 1O 6 M , the 
ionization front quickly expands to a radius of over 1 kpc 
and the halo gas is effectively evacuated with the mean 
density of the halo decreasing to 1 cm~ 3 (see also 
Whalen et al. 2004). In larger mass ( 1O 7 M ) halos, 
the Hn region is confined well inside the virial radius and 
the central density remains high at > 10 4 cm~ 3 . Since 
SN explosions arc triggered at the death of the central 
star, the evolved profiles of density, temperature, and 
velocity of these simulations should serve as appropriate 
initial conditions for the studies of the SN explosions and 
subsequent remnant evolution. 

In the present paper, we study the evolution of SNRs 
in the high-redshift universe, using 1-D numerical sim- 
ulations that start from the final outputs of our earlier 
simulations of the first Hn regions mentioned above. We 
primarily consider point explosions caused by Population 
III SN with a total explosion energy of 10 51 — 10 53 erg. 
We calculate the destruction efficiency and the rate of 
radiative loss from the SNRs to address their cosmolog- 
ical implications. Our study differs from previous works 
on astrophysical explosions in the following points; 1) 
our simulations include the gravitational force exerted 
by the host dark matter halo, and 2) we adopt realistic 
initial conditions for Population III SN sites that are cal- 



culated by self-consistent radiative transfer calculations. 
The simulations incorporate all the relevant cooling pro- 
cesses at T > 10 4 K for a primordial gas, including the 
inverse- Compton cooling which is important at z > 10. 
We also examine several cases with an enormous explo- 
sion energy, E tot = 10 54 erg, on the assumption that mul- 
tiple massive stars form in a single star-forming region 
and that they explode almost simultaneously. This situ- 
ation may indeed be plausible in large mass systems in 
which gas can cool efficiently via hydrogen atomic cool- 
ing. Because the lifetime of massive stars (supernova 
progenitors) is typically comparable to or shorter than 
the evolution time scale of SNRs, we approximate such 
simultaneous multiple explosions by a single explosion 
with a large total energy. 

Throughout the paper, we work with a A-dominated 
cosmology with matter density Om = 0.3, cosmological 
constant Qa = 0.7, the Hubble constant h = 0.7 in units 
of 100 km s _1 Mpc _1 , and baryon density fib = 0.05. 

2. THE SIMULATIONS 
2.1. The code 

We use the hydrodynamics code of Kitayama & Ikeuchi 
(2000) and Kitayama et al. (2000, 2001, 2004). It 
employs the second-order Lagrangian finite-difference 
scheme in spherically symmetric geometry (Bowers & 
Wilson 1991), and treats self-consistently gravitational 
force, hydrodynamics, and radiative cooling. The basic 
equations are presented in Sec. 2 of Kitayama & Ikeuchi 
(2000). We adopt an artificial viscosity formulation of 
Caramana, Shashkov & Whalen (1998), designed to dis- 
tinguish between shock-wave and uniform compression 
using an advection limiter. 

The gas is assumed to have the primordial composi- 
tion with hydrogen and helium mass fractions, X = 0.76 
and Y = 0.24, respectively. The atomic cooling rates 
are taken from the compilation of Fukugita & Kawasaki 
(1994). Since our primary interests lie in the evolution 
of SNRs at temperatures T > 10 4 K until the formation 
of cooled dense shells, formation of molecules and cool- 
ing by them at lower temperatures are not included in 
the present paper to avoid computational difficulties. As 
discussed in Appendix, we can assume that electron-ion 
equilibration is a reasonable assumption in the present 
calculation. The results of basic code tests relevant to 
SNR evolution are presented in Appendix. From the re- 
sults, we conclude that our Lagrangian code accurately 
reproduces analytic solutions and that it also captures 
shock regions remarkably well. 

2.2. Initial conditions 

Supernova explosions are triggered after the progeni- 
tor stars have shined for their lifetime. Consequently the 
ambient gas density and temperature profiles are signifi- 
cantly modified by their radiation. As shown in Paper I, 
the resulting density and temperature profiles are indeed 
quite complex, being dependent on the central source lu- 
minosity, density profile, and the host halo mass. We 
thus take the final output of our previous simulations for 
the ionization front (I-front) propagation as initial con- 
figurations for the present explosion simulations. 

Unless otherwise stated, we adopt the fiducial set up 
in Paper I; a single progenitor star of mass M s = 200M Q 
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Fig. 1. — Expansion of the Hn region around a M B = 200 Mq 
star embedded in a halo with M h = 3.2 X 10 5 M Q and z c = 20; (a) 
hydrogen density, (b) gas temperature, and (c) outward velocity at 
U = (dotted lines), 1.8 X 10 5 (dashed) and 2.2 X 10 6 yr (solid), 
where ti denotes the time elapsed since the birth of the massive 
star. The vertical thin dotted line indicates the halo virial radius. 



is placed at the center of a halo collapsed at redshift 
z c = 20. The gas density profile prior to the birth of the 
star is p oc r~ 2 and is subsequently altered by the stel- 
lar radiation during its lifetime 2.2 Myr (Schaerer 2002). 
Figure n shows a characteristic example of such runs in 
the case of halo mass M h = 3.2 x 10 5 M Q . The whole 
halo is ionized promptly and the gas inside the central 
~ 100 pc is evacuated within the lifetime of the massive 
star, with the inner gas density being ~ 0.1 — 1cm 3 . For 
larger mass halos with Mh £ 10 7 M©, on the other hand, 
the I-front expansion is greatly hindered. Gas infall from 
the outer envelope continues in such cases and the den- 
sity within the compact Hn region remains high at > 10 4 
cm~ 3 (Fig. 2 of Paper I). 

When making the initial conditions from the outputs 
of the Hn region calculations, we need to add extra shells 
and then interpolate the physical quantities in order to 



resolve thin cooling shocks. We adopt the following 
second-order re-zoning scheme so that their monotonicity 
is maintained as well as mass, momentum and thermal 
energy are conserved in each shell. First, gas density 
is expressed as a function of enclosed volume as p(V). 
For the z-th shell, we denote its center of volume by V$ 
and assign the mean gas density of the shell pi to V$, i.e. 
p(Vi) = pi. Second, the density gradient within the ?-th 
shell is defined by 



(dp/dV) l = minmod(s i „ 1 /2,s i+1 / 2 ), 



(1) 



where 



««,i=bK)-«-i)]/W - M-i), (2) 
^^.^WHW**^ (3) 

The gas density at an arbitrary point within the i-th 
shell is computed by p(V) = p{Vi) + {dp/dV) l x (V - Vi). 
Finally, similar procedures are repeated for the internal 
energy per unit mass and for the velocity, except that 
they are expressed as a function of enclosed mass as u(m) 
and v(m), respectively. 

For comparison, we also perform some runs without the 
radiative feedback prior to the SN explosion. Specifically, 
we adopt the same initial configurations as those in Paper 
I; the density profile follows a pure power-law with p cx 
r -2 , the gas is isothermal and infalling with the free-fall 
velocity (dotted lines shown as initial configuration in 
Fig. 2]). We adopt this somewhat artificial configuration 
in order to examine how the final results of our explosion 
simulations are affected by the initial density profiles. 

All the simulations are initiated when the SN ejecta, 
assumed to have M e j = M s , has swept up the surround- 
ing mass equal to its own, i.e., at the end of the initial 
free-expansion phase (see Append ix). This phase ends 
roughly at t sw given in equation (| A5|) after the explo- 
sion. The extent of the inner-most shell is so chosen that 
the enclosed mass is M e j. We then add the ejecta mass 
M e j and the thermal energy E?sn to this cell assuming 
that the ejecta and the swept-up material are fully ther- 
malized. In some cases, the density of the inner-most 
shell is larger than ~ 10 5 cm -3 , its radiative cooling time 
gets smaller than i sw , and the above assumption may not 
be valid (Terlevich et al. 1992; see also Appendix). We 
therefore adopt the following prescription for the initial 
energy input. To the crudest approximation, the energy 
lost before the Sedov phase is E rat i(t sw ) — E(t sw ) x t sw , 
where E(t sw ) is the rate of energy loss in the inner- 
most shell computed under the above assumption. If 
-EgN — -E'rad(isw) is still greater than zero, we add it to 
the inner-most shell instead of Ssn and start the simula- 
tion. If it is less than zero, on the other hand, we regard 
the SN shock as being stalled during the initial phase and 
do not perform the simulation. More detailed discussion 
on the relevant time-scales is presented in Appendix. 

We vary the input explosion energy £?sn over the range 
10 50 — 10 53 erg to account for various possibilities on 
the fate of massive stars; Type II SN, pair-instability 
SN, or by hypernovae (Heger & Woosley 2002; Umeda 
& Nomoto 2002). We also explore the cases where 10 
progenitor stars with M s — 200M Q have produced an 
Hn region and explode with a total energy of Ssn = 10 54 
erg. Since the lifetime of the massive stars, a few million 
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years, is shorter than the characteristic evolution time of 
the SNR, we approximate the multiple SNe as a single 
explosion at the center with i?sN = 10 54 erg. 

The outer boundary is taken at 1 — 100 times the virial 
radius depending on the evolution. Since we do not know 
a priori the extent to which the blastwave can propa- 
gate, this procedure required us to carry out some of the 
simulations more than once. We set progressively larger 
radius for the outer boundary and made sure that it en- 
closes the SN shock after 10 7 yr and also that the bound- 
ary radius does not affect the result within the shock 
radius. Whenever necessary, we extrapolate the gas den- 
sity to the envelope assuming par -2 . We assume that 
the internal energy and velocity are constant in the ex- 
trapolated region in accord with the adopted isothermal 
density profile. The pressure outside the outer boundary 
is set by extrapolating the pressure inside the boundary 
by a quadratic polynomial. 

The gas shells are spaced such that the shell mass 
changes by a constant ratio, less than 1%, between the 
adjacent shells. The total number of shells is N = 1000. 
We have checked that our results remain unchanged when 
the shell number is doubled or halved. 

For the dark matter component, we assume that the 
density profile is given by a Navarro, Frenk & White 
(1997) profile as in Paper I. We follow Bullock et al. 
(2001) to determine the concentration parameter of a 
halo with total mass Mh collapsing at redshift z c , by ex- 
trapolating their formula to the lower halo masses and 
the higher redshifts. The dark matter density is normal- 
ized so that the dark matter mass within the virial radius 
is equal to Mdm = ( 1 — ^b/^m)M\i- As we focus on the 
evolution within 10 yr after the SN explosion, less than 
the halo dynamical time at z — 20, we assume that the 
dark matter density is unchanged throughout the simu- 
lation. 

3. RESULTS 
3.1. Evolution of early SNRs 

As a characteristic example of an early SNR evolution, 
we plot in Figure [21 the radial profiles for the run with 
E SN = 10 51 erg, M h = 3.2 x 10 5 M Q and z c = 20. In 
what follows, we denote by t the time elapsed since the 
end of the initial free-expansion stage. Prior to the SN 
explosion, the central massive star has ionized the gas 
beyond the virial radius and the associated shock has 
swept away the surrounding gas. The central gas density 
is ~ 0.2 cm~ 3 (FigQ at the time when the explosion 
is triggered. This greatly helps the propagation of the 
newly formed SN shock front. It quickly catches up with 
the foregoing shock generated during the Hii region for- 
mation (see the double peaks in density at t = 3 x 10 5 
yr). The halo gas is evacuated almost completely by 
t = 10 7 yr, and the swept-up material resides in a dense 
'shell' at r ~ 300 pc. We found that the evolution of the 
SN blastwave before it merges with the foregoing shock 
is well described by the Sedov self-similar solution for a 
constant ambient density. 

Figure [21 shows the cumulative energy radiated from 
inside the SN shock front and the mass-weighted mean 
temperature and density of the SN shocked region for 
the same run. At t < 10 4 yr, the mean temperature of 
the shocked region is > 10 6 K and the dominant radia- 
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Fig. 2. — Evolution of an early SNR in the case of Esn = 10 51 
erg, M h = 3.2 X 10 5 M Q , and M a = 200 M Q ; (a) hydrogen density, 
(b) gas temperature, and (c) outward velocity at t = (dotted 
lines), 3 X 10 5 (dashed) and 10 7 yr (solid), where t denotes the time 
elapsed since the end of the free expansion stage. The vertical thin 
dotted line indicates the halo virial radius. 



tive processes are free-free and Compton cooling. The 
total emission rate is the la rgest at t ~ 10 5 yr, when the 
dense shell forms (eq. IA12(1 and the mean density of the 
shocked region is the largest. Bulk of the thermal en- 
ergy is then lost via collisional excitation cooling of He + 
and H. Thereafter the SNR enters the pressure-driven 
expansion stage. 

We note here that our code does not solve the radia- 
tive transfer of photons emitted from the remnant itself. 
As the blast wave is induced by the pressure gap much 
greater than that attained by photoionization, the over- 
all effect of the radiative precursor on the dynamics of 
blastwaves is expected to be small. To check this, we 
have performed a test run, constraining the temperature 
of the surrounding gas to be above 10 4 K to approxi- 
mate the pre-ionization ahead of the SN shock. For the 
same set of parameters as in Figure [51 the shock radius 
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line indicates the initial (before radiative feedback) mean density 
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at t = 10 7 yr changes by only ~ 30%. The difference 
will be even smaller for stronger blastwaves with larger 
Esn- Thus we neglect the effect of radiative precursors 
in the rest of the paper. We refer the reader to a detailed 
study of the radiative transfer and the emission spectra 
by Shull & McKee (1979). 

The case of a larger total SN energy _Z?sn = 10 53 erg 
is displayed in Figures 0] and El The gas is quickly evac- 
uated and the shock front exceeds the virial radius at 
t ~ 10 5 yr. The remnant temperature decreases initially 
via adiabatic expansion rather than radiative cooling. 
Due both to high temperature and low density of the 
shocked region, Compton cooling (oc n e T e ) dominates 
over the other radiative processes throughout the evolu- 



FlG. 4. — Same as Figure |2] except for Esn = 10 53 erg. 



tion. In this case, approximately 30% of the explosion en- 
ergy is lost via the Compton cooling. Clearly, explosion 
energy can be efficiently transferred to cosmic microwave 
background (CMB) photons only if the explosion energy 
is large and/or the gas density is small. 

The situation drastically changes if there were no I- 
front expansion prior to the SN explosion. If the ini- 
tial density profile of the run shown in Figure [21 is a 
pure power-law with p cx r~ 2 , the cooling time of the 
inner-most shell gets smaller than t sw and the ejected en- 
ergy is rapidly lost during the free-expansion stage, i.e. 
-E'rad(isw) > ^sn- Accordingly, the blastwave stalls in 
the dense environment and the halo gas in the outer en- 
velope will not be disturbed. This is the case even if the 
ejected SN energy is much greater than the binding en- 
ergy, Sbin ~ 10 erg, for baryons within the virial radius 
of a Mh = 3.2 x 10 5 M Q halo. It clearly shows the im- 
portance of setting-up appropriate initial configurations 
in quantifying the degree of SN feedback. We emphasize 
that simple analytic estimates based on explosion energy 
to binding energy ratio are inappropriate. 
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Corresponding cases of no I-front expansion are real- 
ized in a halo as large as = 10 7 M Q . As shown in 
Figure El the central gas density remains at > 10 5 cm~ 3 
even when the radiation from the massive star has photo- 
ionized and photo- heated the gas (see Paper I). The 
ejected SN energy is lost mainly via free-free emission 
from the shock interior and collisional excitation cooling 
from just behind the shock front (Fig. 0). The latter pro- 
cess efficiently converts kinetic energy into radiation in 
the present case; the gas behind the blastwave is shock- 
heated but radiates quickly at T ~ 10 4 K. This process 
dominates the total cooling rate at t > 30 yr. Due to 
high expansion velocity and high cooling efficiency, an 
extremely dense shell forms at t = 80 yr. The shock ex- 
pands only to ~ 10 pc in t = 2 x 10 6 yr and stalls there- 
after, being well inside the virial radius. Gravitational 
force by dark matter pulls the gas to the center, and the 
mass inflow is eventually recovered. The mass deposi- 
tion rate inside r — 1 pc is ~ 0.1 M Q yr -1 at t = 10 7 
yr. In order to evacuate the gas out of the virial radius, 
multiple SN explosions with _Esn ~ 10 54 erg is required. 



Fig. 6. — Same as Figure 121 except for Esn = 10 53 erg and 
Mh = 10 7 Mq. Note that we plot profiles for a broader range of 
radial distance than in Figure|2]and|l] The peak velocity at t = 80 
yr is 1400 km s — 1 . 



Again this is much greater than the binding energy for 
baryons within the virial radius, E^in ~3x 10 51 erg. 

It is illustrative to compare our 1-D calculation results 
with those of 3-D simulations. We have attempted a 
comparison with the Smoothed Particle Hydrodynamis 
simulations of Bromm et al (2003) as follows. First, we 
set E SN = 10 51 and 10 53 erg in a M h = 10 6 M halo 
collapsing at z = 20 as found in the parent cosmological 
simulation of Bromm et al (2003). Since they used an 
approximate treatment for the radiative feedback prior 
to the explosion, we decided to use the resulting gas den- 
sity profile, rather than directly modeling the radiative 
feedback in their manner. The gas density profile is given 
such that there is a flat core with njj = 1 cm within 
the radius r core — 20 pc and the envelope tih oc r~ 2 be- 
yond that, as described in Bromm et al (2003). Since 
all the gas in the halo is ionized within the lifetime of 
the massive star, we assume that the initial tempera- 
ture is uniform at T — 3 x 10 4 K. To allow sufficient 
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spatial resolution near the center, the mass of the inner- 
most shell is set to be 1 M©, corresponding to the initial 
radius of 2.5 pc for this cell. The shell mass increases 
by 1% between the adjacent shells toward the envelope. 
We ran the simulations starting from this configuration. 
The radial positions of the dense shell are found to be 
r sh = 13, 33, 89, and 214 pc at t = 10 4 , 10 5 , 10 6 , and 
10 7 yr for E SN = 10 51 erg, and r sh = 13, 34, 130 and 570 
pc at t = 10 3 , 10 4 , 10 5 , and 10 6 yr for £ S n = 10 53 erg, 
respectively. Overall the evolution of the shock radius 
agrees well with the results of Bromm et al (2003). At 
r s h £ 30 pc, our calculations predict about a few times 
larger shock radius (i.e. faster expansion) than that plot- 
ted in Fig. 1 of Bromm et al (2003). We have checked 
that our simulations accurately reproduce the analyti- 
cal solutions for the blast wave propagation in the Sedov 
phase (Appendix). Thus the disagreement at the small 
radii could be due to the limited spatial resolution of 
the 3-D simulations, and/or due to details in the way 
the initial explosion energy is deposited. Given that the 



density and velocity fields are likely more complex in the 
3-D simulations, performing more accurate comparison 
would be difficult, and we do not attempt to do so fur- 
ther. We emphasize that it is indeed remarkable that the 
late-time evolution of the SNR is quite similar in both 
cases. 

3.2. Radiative loss from early SNR 

We have shown that the dominant cooling process dur- 
ing the SNR evolution varies with the ambient gas den- 
sity and the explosion energy. In Figure |S1 we plot the 
radiated energy from the SN shocked region within 10 7 
yr as a function of halo mass. The three panels differ 
in the explosion energy, as indicated in the figure. Note 
that, for low-mass halos, the total radiated energy is less 
than E$n because a part of the SN energy is still in the 
form of kinetic energy to drive the shell expansion. We 
find that there are at least three key radiative processes 
in the evolution of high-redshift SNRs. 

First, collisional excitation cooling of hydrogen and he- 
lium is dominant for high Mh or low £?sn (e.g., Figs|3|and 
[7J). It is particularly important in the later stage of the 
SNR evolution, when the temperature drops from 10 7 K 
to 10 5 K and a thin-shell appears. The most prominent 
process is hydrogen Lya emission. For £"sn = 10 51 erg 
and Mhaio = 3.2 x 10 5 M Q , about 25% of the ejected 
energy is ultimately lost by the Lya line in 5 x 10 5 yr, 
yielding an time-averaged luminosity L a ~ 10 erg s~ 4 . 

Second, Compton cooling is dominant for either very 
low Mh or very high -Esn; in this case the remnant ex- 
pands very rapidly and cools mainly after its interior den- 
sity has decreased significantly (e.g., Fig. 0). It is worth 
noting that these SNRs can produce small angular scale 
fluctuations of the CMB via the Sunyaev-Zel'dovich ef- 
fect (e.g. Oh, Cooray & Kamionkowski 2003). In an 
extreme case of .Esn = 10 54 erg and M h = 3.2 x 10 5 M Q , 
~ 30% of the ejected energy is lost by Compton cool- 
ing in 10 7 yr, yielding Eq ~ 10 39 erg s" 1 . The halo 
virial radius corresponds to 0.06 arcsec, and the mean 
y-parameter within it is 3 x 10~ 6 in the first 10 4 yr and 
drops thereafter as the gas temperature decreases. We 
stress that this is a rather extreme case with the ex- 
plosion of ten very massive stars within a small mass 
halo. The total energy lost via Compton cooling declines 
rapidly for larger Mh or smaller £sn- Thus the above 
numbers should serve as the most optimistic estimates 
for the Sunyaev-Zel'dovich effect caused by an individ- 
ual SNRs. 

Finally, free-free emission is important for massive ha- 
los (> 10 7 M o ) in which I-front expansion is suppressed. 
A large fraction of the explosion energy is lost immedi- 
ately in the high density and high temperature region 
near the center, as in the case shown in Figs El and 
[3 Although the mechanical feedback is rather weak in 
this case, these SNRs can be significant compact X-ray 
sources. In the case of Esn = 10 53 erg and Mh = 10 7 M , 
for example, ~ 15% of the ejected energy is lost by free- 
free emission in 30 yr, yielding the rest-frame luminosity 
of Lx ~ 10 43 erg s _1 (Fig. UJ. The corresponding num- 
ber for E SN = 10 51 erg and M h = 3.2 x 10 5 M Q is - 1% 
in 3 x 10 5 yr (Fig. OJl, yielding L x 



10 36 erg s" 1 
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Fig. 8. — The energy radiated from the SN shocked region within 
10 7 yr for Eg-^ = 10 B1 erg (top panel), Eg-^ = 10 53 erg (middle 
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mass. 
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Fig. 9. — -Esn — Afh diagram summarizing the fate of host halos. 
Halos blown away even in the absence of initial I-front expansion 
are marked by triangles, those blown away only in the presence of 
I-front initial expansion by circles, and those not blown away by 
crosses. Dashed and dotted lines show the binding energy of the 
gas for a given and 300 times the same quantity, respectively. 



It is our primary goal in the present paper to quantify 
the halo destruction efficiency by SNe. The fate of host 
halos are summarized in Figure for a wide range of 
halo mass and explosion energy. We regard a halo to 
be 'blown away', if more than 90% of the gas inside the 
virial radius is expelled within 10 yr. We have checked 
that the result is insensitive to the specific choice of this 
threshold between 50% and 99%. 

Figure [5] shows a clear trend that greater SN energy 
is required to blow away larger mass halos. It should be 
noted, however, that the required SN energy is roughly 
300 times larger than the total binding energy for baryons 
in a given halo, if I-front expansion (and hence the de- 
crease of the gas density) prior to the explosion is inef- 
fective. As discussed earlier, this is due to efficient cool- 
ing in the high density (> 10 4 cm~ 3 ) environment near 
the halo center. Once the I-front expansion is properly 
taken into account, low-mass halos with Mh ^ 10 6 Mq 
are blown away even by a single SN with Esn = 10 50 
erg. For higher mass halos (Mh ^ 10 7 M Q ), the impact 
of I-front expansion weakens sharply as shown in Paper 
I. 

We summarize our findings as follows, emphasizing a 
close link between the fate of halo gas classified in Figure 
Eland the main radiative processes described in Sec 13.21 
Compton cooling tends to dominate over other processes 
in halos most readily blown away (triangles). These 
SNRs will cause distortion in the CMB spectrum via the 
Sunyaev-Zel'dovich effect. Collisional excitation cooling 
is likely to be dominant in the others (circles and crosses). 
Ultra-violet photons from these SNRs may contribute to 
reionizing the IGM. Free-free emission is stronger in ha- 
los not blown away by SN explosion (crosses). They may 
be the first X-ray sources in the universe. 

So far we have focused our attention on the evolu- 
tion of SNRs until substantial gas cooling takes place 
within them. It is well known that expanding and cool- 
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ing gas shells are subject to both thermal and gravi- 
tational instabilities (Ostriker & Cowie 1981; Vishniac 
1983; MacLow & Norman 1993). Simple analytic ar- 
guments show that the fragmentation of the swept-up 
shell occurs efficiently when the ambient density is higher 
than ~ 10 cm -3 for the explosion energy of 10 51 — 10 52 
erg (Salvaterra, Ferrara & Schneider 2004). Following 
Machida et al. (2005), we have compared the local free- 
fall time-scale 

(4) 

and the sound crossing time-scale 



where p s h is the mean mass density, Ar s h is the thickness 
and c s is the sound speed, of the thin shell. We define 
the shell in our simulations as the region within the full 
width at half maximum (FWHM) of the density peak, 
and assume that it is stable against fragmentation if tg > 

tsc • 

We find that the above stability condition is satisfied in 
all the 'blown-away' cases with the initial I-front expan- 
sion (triangles and circles in Figure|5J|. In these runs, the 
initial gas density is lower than 1 cm -3 due to radiative 
feedback, resulting in the low efficiencies of both cooling 
and deceleration of the SN shocks. For the run plotted in 
Figure|2(Af h = 3.2 x 10 5 M and E SN = 10 51 erg), the 
shell attains the highest density at t ~ 3 x 10 5 yr with the 
FWHM Ar s h ~ 0.5 pc resolved with eight gas meshes. At 
this point, the above time-scales are <g~2x 10 7 yr and 
t sc ~ 6 x 10 4 yr, respectively. On the other hand, the shell 
may become unstable in some 'stalled' cases, in which the 
shell expansion slows down soon after the internal ther- 
mal energy is lost and cannot effectively counteract con- 
traction. For the run plotted in Figure© (M h = 1O 7 M 
and -Esn = 10 53 erg), tg becomes shorter than t sc when 
an extremely dense shell forms at t = 80 yr. Note that 
including molecular hydrogen cooling would further pro- 
mote fragmention by lowering the temperature and in- 
creasing the density of the shell. Precise account of this 
effect needs further study and we address it in greater de- 
tail in a separate paper (Kitayama & Yoshida, in prepa- 
ration). Finally, it is worth noting that first SN explo- 
sions could trigger self-propagating star-formation (e.g. 
Dopita 1985; Efremov & Elmegreen 1998). In large ha- 
los, second-generation stars could form in SN shells in 
the induced manner, rather than via gas condensation 
onto the halo center. 

4. DISCUSSION 

The first generation stars cause three key feedback ef- 
fects on subsequent star/galaxy formation: (1) Radiative 
feedback by photo-ionizing, -heating, and -dissociating 
the IGM as well as primordial gas clouds, (2) Mechanical 
feedback by expelling the surroundings, and (3) Chemi- 
cal feedback by distributing heavy elements. 

Radiative feedback usually has negative effects on 
star-formation efficiency. The gas in low-mass 
(Mh 10 6 Mq) halos is almost completely ionized even 
by a single massive star (Paper I) . Soft-UV photons in the 
Lyman- Werner bands can dissociate hydrogen molecules, 
the main coolant in the primordial gas, in slightly larger 



mass halos. While Paper I confirmed a positive feedback 
effect on H2 formation just ahead of the I-front (Ricotti 
et al. 2001), it is shown to be rather weak and tentative. 

It has also been argued that soft X-rays increase the 
ionization fraction and promote the formation of hydro- 
gen molecules via the H~ path (Haiman, Abel & Rees 
2001; Oh 2001; Venkatesan, Giroux & Shull 2001). Heat- 
ing by secondary electrons may, however, be substantial 
and the net effect of soft X-rays is still unclear. Indeed, 3- 
D adaptive-mesh-refinement simulations by Machacek et 
al. (2003) showed that, whereas molecular hydrogen for- 
mation is promoted in dense gas clouds that have already 
formed, the dominant effect of the X-ray background is 
to heat the IGM. Machacek et al. also conclude that the 
net effect (either positive or negative) is found to be quite 
small for X-ray intensities below ~ 10~ 23 erg s _1 cm -2 
Hz -1 sr _1 . From our estimates of the X-ray emission 
efficiency of a single SNR, this level of flux is achieved 
if more than a few thousand SN explosions occurred in 
a comoving cubic mega-parsec volume by z ~ 15. For 
reasonable early star-formation models, this rate is im- 
plausibly high, and thus we conclude that global feedback 
effects by X-ray emission from early SNRs are unimpor- 
tant. 

We further argue that the expected intensity of H2 
lines from an early SNR is rather small. Consider as an 
optimistic limit that the ejected energy of i^sN = 10 51 erg 
is used only to heat the ambient gas to T — 8000K, so 
that the atomic cooling is negligible. The heated region 
should enclose the gas mass of 5 x 10 5 M Q or the total 
mass of 3 x 10 6 M Q . It takes ~ 10 7 yr for the shock 
with T = 8000K to propagate over the virial radius ~ 
220 pc of such a halo at z ~ 20. The expected optimal 
luminosity is then Lh2 ~ 3 x 10 36 erg s _1 . In practice, 
some of the gas is heated to higher temperatures, atomic 
cooling will dominate over H2 lines, and the luminosity 
should be lower than this value. 

Mechanical feedback from the first stars is often cited 
as a destructive process in the context of early struc- 
ture formation. The energy released by a single Popu- 
lation III SN explosion can be as large as ~ 10 53 erg, 
which is much larger than the gravitational binding en- 
ergy of a minihalo with mass ~ 1O 6 M . At the first 
sight, this simple argument appears to support the no- 
tion that high-z SNe are enormously destructive. As we 
have explicitly shown, however, the evolution of the cool- 
ing remnants crucially depends on the properties of the 
surrounding medium, particularly on the central density 
after the gas is re-distributed by radiation. 

The efficiency of radiative feedback declines sharply 
above a certain halo mass Mh ~ 10 7 M Q (Paper I), 
so does that of the mechanical feedback by SNe. Ha- 
los above this threshold can thus be ideal candidates 
for the hosts of the first massive black-holes in cases 
where supernovae leave remnant black-holes. Our calcu- 
lations indicate that the mass deposition rate at parsec 
scales reaches ~ O.lM yr _1 in the case of i?sN = 10 53 
erg, M h = 10 7 M Q at t = 10 7 yr (Sec. E3J). This 
value is much higher than a simple estimation for the 
Bondi accretion rate 4 x 10~ 8 (M B h/100M o ) 2 M yr -1 
for the ambient gas density 10 4 cm~ 3 and temperature 
10 4 K, or the rate resulting in the Eddinton luminosity, 
2 x lO- 6 (e/O.l)- 1 (M BH /lOOM )M yr" 1 , where M BH 
is the black-hole mass and e is the efficiency factor of 
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fueling. Therefore, it should at least provide a necessary 
(obviously not sufficient) condition for the near-Eddinton 
growth of black-holes at the center of these large halos 
even after the mechanical feedback from the first SNe. 
We note, however, that actual mass accretion around 
black holes, occurring at much smaller radius than our 
simulations probe, is much more complicated because of 
radiation- hydrodynamical effects, non-spherical geome- 
try, etc. These issues are still highly uncertain and be- 
yond the scope of the present paper. 

As has been often discussed in the literature (Ricotti 
& Ostriker 2004; Yoshida et al. 2004), the first star 
formation is likely to be episodic, at least locally, be- 
cause the feedback effects tend to quench further star- 
formation in the same place. Although metal-enrichment 
by the first stars could greatly enhance the gas cooling 
efficiency, which would then change the mode of star- 
formation to that dominated by low-mass stars (Mackey, 
Bromm & Hcrnquist 2003; Bromm & Loeb 2003), the 
onset of this 'second-generation' stars may be delayed 
particularly in low-mass halos. Our results imply that 
the strength of this negative feedback substantially de- 
creases for £ 10 7 M Q halos. In CDM models, halos grow 
hierarchically and thus all the halos with £ 10 7 M Q used 
to be 'mini-halo' at some earlier epoch. If star-formation 
is prevented earlier in the progenitor small halos either by 
rapid mass accretion (Yoshida et al. 2003) or by a soft- 
UV background which raises the minimum mass scale for 
the gas cloud formation (Machacek et al 2001; Wise & 
Abel 2005), primordial stars could form in £ 10 7 M Q ha- 
los. Only if multiple explosions take place after a burst of 
star- formation in the large halos, a complete destruction 
can take place. Wada & Venkatesan (2003) showed that 
such multiple explosions can easily disrupt pre-galactic 
disks. 

5. CONCLUSIONS 

First supernova explosions in the early universe have 
a number of consequences on subsequent star formation. 
While they tend to expel gas from low-mass halos, the 
extent to which such feedback can operate depends sen- 
sitively on the initial configuration of the medium sur- 
rounding the explosion site. We have shown explicitly 
that the ionization front expansion around a massive 
progenitor star can significantly aid gas evacuation by 
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supernova blastwaves. For halos with mass 1O 6 M0, 
essentially all the gas is expelled by merely a tenth of 
nominal SN explosion energy. On the other hand, if the 
Hll region is confined within the central portion of the 
halo, the gas density remains high, and the necessary SN 
energy for 'blow-out' is larger by more than two orders 
of magnitude than the total binding energy for baryons 
in the halo, i.e., E SN > 10 54 erg, for M h '= 1O 7 M at 
z = 20. This is due to very high efficiency of radiative 
cooling in the dense gas clouds. 

We have also found that there are mainly three radia- 
tive processes that are important in the early SNRs. (i) 
Collisional excitation cooling, particularly via the Lya 
line emission, is dominant for either higher-mass halos 
or lower SN energy, (ii) Compton cooling is dominant 
for either lower-mass halos or higher SN energy, (iii) 
Free-free emission gets stronger for massive (> 10 7 M Q ) 
halos, in which the expansion of an Hn region is greatly 
suppressed, and can make them compact X-ray sources. 

Overall, the emergence of the first generation stars 
have significant impacts on the thermal state of the IGM 
in the early universe in at least two ways; the initially 
neutral cosmic gas is first photo-ionized and heated to 
~ 10 4 K by stellar radiation, and subsequently heated 
to even higher temperature and enriched with heavy ele- 
ments by SNe. Our results imply that radiative and SN 
feedback processes are closely related to each other and 
thus need be treated in a self-consistent manner. Be- 
cause of the strong feedback effects, early star-formation 
is likely to be self-regulating; if the first stars are massive, 
only one period of star-formation is possible for a small 
halo and its descendant. The sharp decline in the effi- 
ciency of both feedback effects at > 10 7 M Q indicates 
that the emergence of large halos can drastically raise the 
global star formation rate of the Universe. This may lead 
to accelerating the metal enrichment of the IGM at an 
early stage of galaxy formation. 
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APPENDIX 

BASIC FEATURES OF THE SNR EVOLUTION IN THE EARLY UNIVERSE 

The evolution of a SNR may be divided into three stages: (1) a free-expansion phase when the SN ejecta expand 
at nearly a constant velocity, (2) an adiabatic Sedov phase which is described by a self-similar solution, and (3) a 
pressure-driven snow-plow phase when a dense shell forms within which the thermal energy of the gas is radiated away. 
In the following, we describe basic elements of physical processes associated with point explosions. 

The free-expansion phase lasts until the ejecta sweeps up roughly the same amount of mass as their own, M C j, in 
the surrounding medium: 
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where R sw is the radius enclosing M e j, m p is the proton mass, X is the hydrogen mass fraction, and «h is the mean 
hydrogen number density of the surrounding medium. With the ejected kinetic energy Esn, the initial shock velocity 
is 
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corresponding to the initial shock temperature of 
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where /i is the mean molecular weight and k is the Boltzmann constant. The time taken for the shock to reach R sw is 
then estimated as 
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The time required for electrons to attain the shock temperature via Coulomb collisions is roughly (Masai 1994) 
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So long as t s < t sw , one can assume that electron-ion equilibration is achieved and T c ~ T s j, where T c is the electron 
temperature. The hot gas in the interior of high redshift (z > 10) SNRs cools initially via free-free emission and 
Compton scattering with the time-scales given respectively by 
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If either tg or ^comp is shorter than i sw , the ejected energy is radiated away before entering the Sedov phase. If ts 
and tcomp are both greater than t sw , the evolution is adiabatic. Note that t s can be larger than t sw for low nn or high 
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Fig. A10. — Propagation of the blastwave with Esn = 10 51 erg in a uniform medium corresponding to a virializcd halo with = 10 6 
Mq and z c = 20. Dashed lines indicate the simulation results at 10 4 and 10 6 yr, while dotted lines denote the corresponding self-similar 
solutions. 



Esn- The radiative cooling, however, is unimportant in such cases and the shock can still be well described by the 
adiabatic solution and electron-ion equilibration is eventually achieved. 

For a power-law density profile with riu oc r~ w , the expansion of a blastwave in the adiabatic Sedov phase is described 
by the self-similar solution, 

/ tp j-2 \ 1/(5-™) 



R s (t) cx , (A9) 

V rin J 

where Rs(t) is the shock radius at time t after the explosion (Sedov 1959; see also Ostriker & McKee 1988). In reality, 
the gas density profile can be more complex because photo-evaporation due to radiation from the central star is very 
effective in low- mass halos (Paper I). If the size of the Hii region is comparable to the virial radius (e.g. Fig. [2Jl, the 
gas density inside the Hii region is nearly uniform at hh ~ 0.3[(1 + z)/20] 3 cm~ 3 . For a constant density medium, the 
Sedov- Taylor self-similar solution reduces to 

{^kT (ik) 2 ^ <™> 

^K & )-( T |f^) 2/5 ( T J-)- 8/5 , ( An, 

Note that the shock radius in equation (|A9() scales only weakly with the density and the explosion energy unless the 
density profile is very steep. The 3- D simulations of Bromm et al. (2003) show that the evolution of the shock radius 
is well described by equation (|A10|) in a piecewise manner for regions having density profiles with different slopes. 
The above adiabatic solution breaks down once the gas becomes radiative. A dense shell will form roughly at the 
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Fig. All. — Same as Figure lAlOl except that the medium has a gas density gradient with index w = 2. 



sag time (Cox 1972; Shull 1980) at which t ~ i ff /10 with T e = T s (t) from equations (|X7|) and (|HT|l : 

The shock will then follow the snow-plow solution, R s (t) oc t 2 / 7 (Shull 1980). The momentum-driven slow expansion 
is still important over a cosmological time-scale (Bertschinger 1983), with the radius doubling every twenty folds in 
time. 



CODE TEST 

We have tested and verified the validity of our Lagrangian code by a variety of test problems. The results of the 
'Stromgren sphere' test is presented in Paper I. We here describe the results of the Sedov explosion problem (Sedov 
1959), which is the most relevant to the present paper. 

We solve the evolution of a blastwave in a spherical halo with Mh = 10 6 M Q and z c = 20. For the sake of direct 
comparisons with the self-similar solutions, we neglect gravitational force and radiative cooling. The initial gas density 
profile is a power-law uh oc r~ w with either w = (uniform) or w — 2. The total number of gas shells is N = 1000; 
the inner-most shell encloses the gas mass 1 M© and the shell mass increases by a constant ratio, 1%, between the 
adjacent shells toward the envelope. The thermal energy of -Esn = 10 51 erg is given to the inner-most shell. The initial 
thermal energy in the other shells is essentially zero. The initial velocity is also taken to be zero. The simulation is 
otherwise perfo rm ed in the same way as the other runs in the present paper. 

Figures TA 1 01 and | A 1 1 1 compare the outputs of the above test runs and analytical solutions (Ostriker & McKee 1988). 
The numerical results reproduce well the latter for both w = and w = 2. The positions of shock front, defined as 
the peak of thermal pressure, agree with the analytical expectations within 1%. This further ensures the capability of 
the code for the simulations in the present paper. 



